home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / ranlib / genmn.f < prev    next >
Text File  |  1996-07-19  |  2KB  |  83 lines
  1.       SUBROUTINE genmn(parm,x,work)
  2. C**********************************************************************
  3. C
  4. C     SUBROUTINE GENMN(PARM,X,WORK)
  5. C              GENerate Multivariate Normal random deviate
  6. C
  7. C
  8. C                              Arguments
  9. C
  10. C
  11. C     PARM --> Parameters needed to generate multivariate normal
  12. C               deviates (MEANV and Cholesky decomposition of
  13. C               COVM). Set by a previous call to SETGMN.
  14. C               1 : 1                - size of deviate, P
  15. C               2 : P + 1            - mean vector
  16. C               P+2 : P*(P+3)/2 + 1  - upper half of cholesky
  17. C                                       decomposition of cov matrix
  18. C                                             REAL PARM(*)
  19. C
  20. C     X    <-- Vector deviate generated.
  21. C                                             REAL X(P)
  22. C
  23. C     WORK <--> Scratch array
  24. C                                             REAL WORK(P)
  25. C
  26. C
  27. C                              Method
  28. C
  29. C
  30. C     1) Generate P independent standard normal deviates - Ei ~ N(0,1)
  31. C
  32. C     2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM
  33. C
  34. C     3) trans(A)E + MEANV ~ N(MEANV,COVM)
  35. C
  36. C**********************************************************************
  37. C     .. Array Arguments ..
  38.       REAL parm(*),work(*),x(*)
  39. C     ..
  40. C     .. Local Scalars ..
  41.       REAL ae
  42.       INTEGER i,icount,j,p
  43. C     ..
  44. C     .. External Functions ..
  45.       REAL snorm
  46.       EXTERNAL snorm
  47. C     ..
  48. C     .. Intrinsic Functions ..
  49.       INTRINSIC int
  50. C     ..
  51. C     .. Executable Statements ..
  52.       p = int(parm(1))
  53. C
  54. C     Generate P independent normal deviates - WORK ~ N(0,1)
  55. C
  56.       DO 10,i = 1,p
  57.           work(i) = snorm()
  58.    10 CONTINUE
  59.       DO 30,i = 1,p
  60. C
  61. C     PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky
  62. C      decomposition of the desired covariance matrix.
  63. C          trans(A)(1,1) = PARM(P+2)
  64. C          trans(A)(2,1) = PARM(P+3)
  65. C          trans(A)(2,2) = PARM(P+2+P)
  66. C          trans(A)(3,1) = PARM(P+4)
  67. C          trans(A)(3,2) = PARM(P+3+P)
  68. C          trans(A)(3,3) = PARM(P+2-1+2P)  ...
  69. C
  70. C     trans(A)*WORK + MEANV ~ N(MEANV,COVM)
  71. C
  72.           icount = 0
  73.           ae = 0.0
  74.           DO 20,j = 1,i
  75.               icount = icount + j - 1
  76.               ae = ae + parm(i+ (j-1)*p-icount+p+1)*work(j)
  77.    20     CONTINUE
  78.           x(i) = ae + parm(i+1)
  79.    30 CONTINUE
  80.       RETURN
  81. C
  82.       END
  83.